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Abstract 

We use the recently introduced nested sampling Monte Carlo algorithm to calculate the 
partition function of the three dimensional periodic hard sphere system. Given the partition 
function, we evaluate the compressibility and the free energy as functions of the packing frac- 
tion and local order, showing that the transition to crystallinity is barrierless, and that the 
entropic contribution of the jammed states to the free energy is negligible for packing fractions 
above the phase transition. We also identify quantitatively the extent of the region of jammed 
states, and find that the maximally random jammed configuration is surprisingly disordered. 

Introduction 

One of the first material systems ever studied using computer simulation was the hard sphere 
model and its simpler two-dimensional version of hard disks J-^ Since then hard spheres have been 
widely used as model particles not just in the field of statistical mechanics, but also in biology, 
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engineering and material science, and consequently extensively studied by analytical means, by 
computer simulations, and also by experimental methods using colloids. 

The hard sphere system undergoes a first order solid-fluid phase transition. 4 The fluid phase is 
the stable phase from packing fraction = up to the freezing point, w 0.494, beyond which 
fluid and solid coexist until the melting point of w 0.545 is reached. 4 5 Above this packing frac- 
tion only the solid phase is stable, and the densest possible packing fraction which can be produced 
with identical hard spheres is n/3\/2 ps 0.74048, corresponding to the close packed face-centered 
cubic (fee) lattice or its stacking variants. The last assertion was known as Kepler's conjecture and 
was recently proved. ^However, metastable solid structures can be produced (both experimentally 
and in simulations) where the particles are disordered at relatively high density, in the so-called 
random close packed (RCP) state. Due to the proposed similarity of this disordered phase to 
glassy phases in other systems, these structures are widely studied, however there is no true glass 
transition in one-component monodisperse hard sphere systems It is known from experimental 
and computer simulations as wellP^ that the disordered phase is thermodynamically 
unstable with respect to the ordered crystal. 

Although the term "random close packed" is widely used in the literature, Torquato and his 
co-workers pointed ouP that the RCP state is very ill defined, and the corresponding packing 
fraction that is reported varies in the range of ~ 0.60 — 0.68^^ depending on the protocol used 
to produce the configurations. The problem lies in the fact that it is not well defined how random 
the 'random', or how close packed the 'close packed' really is. They suggest that the randomness 
should be quantified by an order parameter, and the term 'close packed' should be replaced by the 
term 'jammed' . A system is jammed, if all the particles are jammed in the system, and a particle is 
jammed if it cannot be translated while the rest of the system is kept fixed. The various states can 
be sketched on a two dimensional diagram using the packing fraction and an (as yet undefined, 
"ideal") order parameter *F as axes, as shown in Figure [TJ reproduced form Torquato et al. . The 
bottom left corner (0 = 0, \F = 0) represents the completely disordered spheres of zero size, while 
the top right corner (0 ps 0.74048, *P = 1) corresponds to the perfectly ordered and close packed 
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Figure 1 : Schematic plot of an arbitrary order parameter, *F, of hard sphere configurations as a 
function of the packing fraction, 0. The grey area represents the inaccessible region, and the 
black squares show the maximally random jammed state (MRJ), the jammed state with the lowest 
possible packing fraction (A) and the perfect crystalline structure (B). 

fee crystal. The grey area represents the inaccessible region, i.e. where the hard spheres would 
overlap and all the jammed configurations are located within the loop. Two specific points of the 
loop boundary can be identified: the point marked with A corresponds to the jammed structure 
with the lowest possible packing fraction, and the maximally random jammed (MRJ) state is the 
most randomly packed jammed structure (i.e. having the lowest order parameter). 

In this paper, we report on simulations which sampled the entire phase space of the finite size 
hard sphere system with periodic boundary conditions. This enables us to quantitatively substanti- 
ate the sketch of Figure [1} estimate the parameters of the extremal points of the jammed phase and 
also calculate the free energies of the various states as a function of the packing fraction. 

Thermodynamics of hard spheres 

In the next section we will give a summary of our sampling method that allows the direct compu- 
tation of the partition function. Because this requires us to express the partition function in a form 
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that is amenable to this highly novel approach, we briefly recall the fundamental thermodynamics 
of hard spheres. The isotherm-isobar partition function, A, is given by 



A(/3,p,^) = p^/3p|exp[-/3(//(p,q,y)+py]JpJqjy, (1) 

where p and /3 are the pressure and the inverse thermodynamic temperature, respectively, V is the 
volume, N is the number of particles, H is the Hamiltonian, and p and q are the momentum and 
spatial coordinates. Integrating out the momentum dependent part, Z p , in the usual way, we have 

^ - UwY ,2 < 

A(P,p,N) = PpZ P iP) fdV f dqexp[-P(U(q,V) + pV\. (3) 

Jo Jv N 

We now change to fractional coordinates s: 

A(P,p,N)=l3pZ p (P) [ ds rdVV N exp[-P(U(s,V)+pV]. (4) 

J(0,l) 3N Jo 

For hard spheres of a given radius, the potential energy is either zero or infinity that, for a given set 
of fractional coordinates, depends on the volume, 

ifV>V c (s) 

U(s,V) = { ~ (5) 

°o ifV<V c .(s) 

where V c (s) is the function that gives the lowest permissible volume for the configuration repre- 
sented by the fractional coordinates s. Substituting into eq. 4, 

A(/3,p,iV)=/3pZ p (/3) / ds T dVV N &Lp[-p P V], (6) 

J(0,l)™ JV c (s) 
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where the last integral can be transformed to a standard exponential integral, as 

f'OO 

/ dV V N cxp[-l3pV} = V c (s) N+l E_ N (pv c (s)p). (7) 



poo g Xt 

E N {x) = J dt-p- (8) 
The partition function thus becomes 

A( J 8, J p,/V)= J 8pZ p ( J 8) / dsV c ( S ) N+l E_ N (pv c (s)p) (9) 
= J^ Z P^\f {0l)3 / s nN+l,ms)p), (10) 

where T is the upper incomplete gamma function, 

POO 

T(N,x) = / t N - l e- t dt = x N Ei_ N (x). (11) 

J X 

For completeness, recall that given the partition function, thermodynamic observables can be 
computed easly. For example, the Gibbs function is given by 

G{T,p,N) = -kT\nA, (12) 

and so the expected volume and the compressibility are, 

v = f = Mf < 13 > 

dp A dp 

1 dV 

K = ~vw (14) 
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Nested Sampling of hard spheres 

Nested sampling (NS) was recently introduced by SkillingpSED m the field of applied probabil- 
ity and inference, to sample probability densities in high dimensional spaces where the regions 
contributing most of the probability-mass are exponentially localised. In essence, nested sam- 
pling generates a set of samples (using an algorithm outlined below) that are suitable for the direct 
evaluation of the density of states and thus the partition function at a wide range of its thermody- 
namic parameters and all other observables are then calculated using this discrete approximation. 
A unique strength of the method is that the ampling efficiency is not curtailed by the presence 
of phase transitions. In our previous work^£3| we g ave a detailed description of the nested sam- 
pling method and used it to study the phase space of clusters of Lennard- Jones particles, showing 
enormous gains in efficiency in computing heat capacity curves. Due to its ability to estimate the 
partition function, nested sampling permitted the study of phase stability, and both relative and 
absolute free energies. Since then, Wheatley and coworkers used a sampling scheme that is essen- 
tially identical to nested sampling to study vapour- liquid coexistence in a variety of fluids. EH3 

For the present case of hard spheres, we are going to approximate A numerically by discretising 
the positional integral in eq. 10, 

a* = £w,r(^+ 1,^(80), (15) 

where i runs over the samples, and vv ; is the weight of each sample, with = 1. The algorithm 
generates a series of configurations with decreasing volume levels V c (si) and corresponding 
weight factors w; which are also decreasing: configurations with smaller volumes (i.e. larger 
packing fractions) represent smaller phase space volumes and thus a smaller contribution to the 
partition function. 

1. Start by initialising a set of K configurations by drawing samples from a uniform distribution 
in a very large volume (so that they are all valid configurations), and determine the critical 
volumes V c for all samples. 
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2. Choose the sample point with the largest critical volume V c (i.e. smallest packing fraction). 
In the z'th iteration of the loop, this gives the V c (si) for the sum in eq. 15. (We explain below 
how to compute w ( -.) 

3. Remove this sample point and replace it by duplicating another one which is randomly cho- 
sen from the rest of the existing sample set. 

4. Perform a random walk on the newly created configuration, with both the positions of the 
particles and the shape of the simulation box changing, but the critical volume V c is never 
allowed to exceed that of the configuration selected in step 2. The walk is performed for 
a fixed number of steps, M, the aim being to decorrelate the newly created configuration 
from its progenitor, and thus end up with K samples that are now uniformly distributed with 
critical volume V c (s) < V C (S(). 

5. Go to step 2. 

Because in each iteration we throw away the single sample with the highest critical volume 
from a set of K that are uniformly distributed below a given critical volume, the remaining sample 
set represents a phase space volume that is a factor [K—\)/K smaller than that of the previous 
set, so the phase space volume corresponding to the sample set at the z'th iteration is [(K— l)/K] 1 . 
Thus, a sample point that was removed in the z'th iteration contributes to the overall sum with a 
weight wi that is the difference between the phase space volumes corresponding to the two sets 
before and after it was thrown away, w ( - = [(K - 1 )/£]'' - [(K - \)/K] i+l . Since the sampled 
phase space volumes decrease exponentially with iteration number, the overall sampling converges 
quickly, and is able to locate exponentially small parts of configurational space. At the same time, 
because in each iteration the sampled phase space volume decreases by a factor of (K—l)/K 
that is fixed (even near V c values that correspond to discontinuous phase transitions), there is no 
inherent difficulty in generating the samples, there is no "equilibration problem": if the phase space 
volumes decreases rapidly with V c , the successive V c (Si) values will be closely spaced. A larger 
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K leads to slower convergence but correspondingly higher resolution in the discretisation of the 
partition function. 

Also note the absence of the temperature /3 or pressure p in the actual sampling algorithm. 
Since T(N + l,j3pV c ) that appears in eq. 15 is a monotonic function of V c , the above derivation 
of the sampling weights is independent of /3 p. Thus the expectation value of any observable can 
be evaluated at an arbitrary pressure just by resumming over the generated samples, obviating the 
need to perform a separate simulation for each desired pressure. 

The size of sample set, K, and the number of attempted moves, M, during the random walk in 
step 4 of the algorithm are convergence parameters. In the present work, K was between 2000 and 
5000, and M was between 160 and 480. In order to quantify finite size effects, systems of 32, 64 
and 128 hard spheres were used. For each system size we performed four independent calculations 
starting from different initial conditions. The total number of MC steps needed to converge the 
partition function and to reliably reach the crystalline structure with the highest possible packing 
fraction are shown in Table [T] In comparison, Noya et. aP performed 4 x 10 9 — 5 x 10 10 MC 
steps to converge a coexistence simulation at a single packing fraction, using a unit cell of 5184 
particles. 

Table 1 : Number of MC steps needed to produce a converged nested sampling run and to reliably 
reach the crystalline structure with the highest possible packing fraction in case of different system 
sizes. 



Number of particles 


Number of MC steps 


32 


9.0 x 10 y 


64 


2.0 x 10 10 


128 


7.1 x 10 10 



Compressibility 

Given the discrete approximation of the partition function, it is easy to evaluate thermodynamic 
observables, such as the compressibility as a function of packing fraction, which is shown in Fig- 
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Figure 2: Compressibility (in units of /3a 3 ) as a function of packing fraction for 32 (red), 64 
(green), 128 (blue) and 128 with longer random walks (magenta) hard sphere particles in a periodic 
simulation box. The standard errors are shown by dashed lines. 

ure[2j There is a well defined peak at = 0.49 — 0.55, which matches the literature value for the 
phase transition region.^ The compressibility clearly has a significant finite size error for N = 32 
particles, but it is remarkable that the finite size error of the nested sampling simulation is very 
small already for N = 64. 

Free energy 

Access to the partition function enables the calculation of absolute free energies as a post-processing 
step after the nested sampling simulation is complete, therefore we are able to compute the Gibbs 
free energy of the hard sphere system. The natural variable of the Gibbs free energy is the pres- 
sure, and to shed light on the phase transition, we introduce the well known bond order parameter 
Qj^as a further thermodynamic variable. Its value is zero for completely disordered systems, and 
increases with degrees of crystallization, however the Qd value differs for the perfect fee and hep 
(and their stacking variants) crystals, Qg(fcc) = 0.57452 and <2 6 (hcp) = 0.48476. Thus the free 
energy becomes 
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G(p, fie) = ln(A(j8,/>, fie)) « -/3 1 In 



£(w ; r(^V+ l,j8/»V c (si))fi(fi 6 (si) - fie)) 



, (16) 



where 8 is the Kronecker delta function. The calculated Gibbs free energies are shown in Figure [3] 
as a function of Q(, for a series of fixed pressure values (with the corresponding average packing 
fraction in parentheses). The phase transition at p m 1 1.2, rs 0.5 is seen to be almost barrierless. 
Above the phase transition, two distinct free energy minima appear, corresponding to the hep and 
fee stacking orders, with a barrier between them. 



Order parameter vs packing fraction 

To quantitatively substantiate the diagram of Figure [T] we again use the bond order parameter fig 
as the order parameter for the vertical axis. The calculated fig parameter of the nested sampling 
configurations as a function of their packing fraction are shown with black dots in Figure [4] for 
all system sizes. The phase transition corresponding to the compressibility peak near <j> w 0.49 is 
revealed as a sharp increase in crystallinity. 

Nested sampling allocates samples to various regions of phase space in proportion to the re- 
gion's phase space volume. Thus white areas in Figure [4] (i.e. regions with no samples in our 
simulation) correspond to regions with very small phase space volume (approximately < 1/K for 
a given packing fraction). Accordingly, at low packing fraction, the crystalline states while clearly 
accessible, occupy a negligible phase space volume. Conversely, at high packing fraction, the 
disordered states have negligible volume, which manifests itself in the well known result that the 
disordered phases are thermodynamically unstable here. The samples with the highest possible 
packing fraction fall into distinct groups with different order parameter values, because the nested 
sampling simulation found both the fee and the hep phases. However, it is interesting to see that 
the branches leading to either stacking variants of the crystal separate at relatively small <j) values, 
already around the melting point. 
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Figure 3: Gibbs free energy per particle of a periodic lattice of hard spheres with 128 particles 
in the unit cell as a function of the order parameter at different pressure values (/3p in o 
units) calculated by eq. 16. In order to show the curves for different pressure values on the same 
plot, p(V) was added to the free energy. The corresponding packing fraction (0) values are also 
indicated. The bottom three panels show the free energy curves in a larger scale around the phase 
transition. 
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Figure 4: Order parameter as a function of the packing fraction for Af = 32,64 and 128 hard 
spheres. The black dots represent the configurations generated during the nested sampling runs 
and the coloured symbols represent the quenched jammed configurations (see text). The last panel 
shows the area occupied by the jammed states of 32 (red), 64 (green), 128 (blue) using a contour 
line after gaussian smearing. Perfect fee and hep structures are indicated by filled black circles, 
and the points MRJ (maximally random jammed state) and A (jammed state with lowest possible 
packing fraction) are also shown on the last panel for the 128 particle system. 
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To distinguish the region of jammed states from the totally inaccessible region of high packing 
fraction and low order at the right hand side of the plot, we chose 10000 points from the nested 
sampling configurations in the packing fraction range of 0.40 — 0.73 and rapidly quenched them 
by increasing the particle diameter (and hence the packing fraction) with minimal change in the 
location of particle centers until the system became jammed. 

The final jammed configurations are plotted by coloured stars in Figure |4j The size and shape 
of the region they occupy is converges with increasing particle number to the region of RCP states. 
It can be seen from these figures that jammed configurations exist in states corresponding to a wide 
variety of order and packing fraction. There are jammed states with packing fraction smaller than 
0.6, and many others are as disordered as the fluid phase, or almost perfect crystals with different 
defects. 

The quenching process allows the characterisation of the phase transition without the use of 
an explicit and somewhat arbitrarily chosen order parameter. In Figure [5] we plot a random subset 
of our samples with the packing fraction before and after quenching on the horizontal and vertical 
axes. The phase transition is again visible, showing that configurations just below and above the 
critical packing fraction are close to jammed states with relatively low and high packing fraction, 
respectively. 

Conclusion 

We have studied the phase transition of the periodic hard sphere system in detail using nested 
sampling. A remarkable property of nested sampling is that the finite size effect is so small that the 
compressibility as a function of packing fraction converges very well over the entire range already 
for 128 particles. Coexistence simulations using the same number of MC moves and ~ 50 times 
more particles obtained just the location of the transition.^ 

Using the samples collected, we were able to plot the free energy as a function of order param- 
eter and packing fraction, showing the phase transition explicitly and quantitatively, demonstrating 
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Figure 5: Packing fraction of the quenched jammed configurations as a function of their initial 
packing fraction for the 32, 64 and 128 hard spheres systems. 

that the transition to crystallinity is barrierless, and that above the phase transition the entropic 
contribution of the jammed states to the free energy is negligible. The global unbiased nature of 
this sampling scheme allowed us to study the phase diagram using the packing fraction and the 
Qd order parameter as thermodynamic variables, identifying the extent of the jammed phase. The 
"maximally random jammed" state, i.e. the jammed state with the lowest £>6 order parameter is 
surprisingly disordered, with a value of less than 0.05. 
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Supplementary 

During the quenching process the size of the diameter was increased at every step until the whole 
configuration was jammed. We considered a configuration jammed when less then 10% of the 
particles could be moved by a maximum length of 0.01% of the sphere radius, thus the trial move- 
ments in the rest of the system were rejected. This means that the packing fraction was maximised, 
further relaxation could only increase it by less then 0.02% in 60000 trial movements of the parti- 
cles. 

The question arises how the location of a given jammed state depends on its initial configuration 
it was minimised from. We chose several packing fraction values where we picked the current 
sample set of the nested sampling run in order to have sets of points uniformly distributed at given 
packing fraction levels and we quenched those with the same method previously described. The 
chosen packing fraction levels and the quenched jammed states are shown with matching colours 
in Figure [6] for a 64 and two separate runs of the 128 hard spheres systems. The figures show the 
same general picture for both system sizes. Different initial packing fraction levels result different 
groups of jammed states, starting from the fluid phase (<j) ~ 0.45) the jammed states have relatively 
small packing fraction and low order, while in case of quenching configurations above the freezing 
point the jammed states end up in the region of crystalline structures. 

This suggests that around the freezing point the relative phase space volume of the crystalline 
structures is small, approximately only 0.2% of the points quenched from the set w 0.485 ended 
up with a high order parameter. Thus, a large part of the phase space here leads to different 
disordered local minima. However, above the freezing point, the relative phase space volume 
of ordered structures becomes dominant, no matter that an immense number of distinct basins 
(corresponding to jammed states with disordered structures) exist at the same packing fraction 
level. 
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Figure 6: Order parameter as a function of the packing fraction in case of a 64 hard sphere sys- 
tem (top panels) and two of the 128 hard spheres systems (bottom panels). The red dots represent 
the configurations generated during the nested sampling runs. The vertical lines indicate the pack- 
ing fractions at which points the sample sets were saved and later quenched in order to get jammed 
configurations. The quenched jammed configurations are shown by symbols of matching colours. 
In the top right panel some of the initial and their corresponding jammed structures are connected 
by green lines to show how the quenching changes the order parameter. The points corresponding 
to perfect fee and hep structures are shown by black circles. 
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